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Abstract. Recent observations of NGC 1068 and other AGN support the idea of a geometrically and optically thick dust torus 
surrounding the central supermassive black hole and accretion disk of AGN. In type 2 AGN, the torus is seen roughly edge-on, 
leading to obscuration of the central radiation source and a silicate absorption feature near lOyt/m. While most of the current 
torus models distribute the dust smoothly, there is growing evidence that the dust must be arranged in clouds. We describe a 
new method for modeling near- and mid-infrared emission of 3-dimensional clumpy tori using Monte Carlo simulations. We 
calculate the radiation fields of individual clouds at various distances from the AGN and distribute these clouds within the torus 
region. The properties of the individual clouds and their distribution within the torus are determined from a theoretical approach 
of self-gravitating clouds close to the shear limit in a gravitational potential. We demonstrate that dumpiness in AGN tori can 
overcome the problem of over-pronounced silicate features. Finally, we present model calculations for the prototypical Seyfert 
2 galaxy NGC 1068 and compare them to recent high-resolution measurements. Our model is able to reproduce both the SED 
and the interferometric observations of NGC 1068 in the near- and mid-infrared. 



Key words. Galaxies: Seyfert - 
- Radiative Transfer 



1. Introduction 



The Unification Scheme for AGN was introduced by 
iMiller & Antonuccil (ll983) to explain the different appearances 
of various types of AGN. According to our present understand- 
ing, the centers of all AGN host a central black hole of several 
10^ to 10^^ Mq. It is surrounded by an accretion disk which 
reaches inwards to the closest stable Keplerian orbit. Due to 
energy gained from accretion in a gravitational potential, the 
accretion disk is heated up to -10^ K, which results in strong 
emission in the UV/optical part of the spectrum. Such a UV- 
bump has been observed in QSOs as well as in Seyfert 1 galax- 
ies. 

At distances where the ambient temperature falls below 
~ 10^ K, dust is able to survive; it is thus possible to have a large 
dusty structure present which surrounds the AGN. This dust 
reservoir can explain the optical properties of type 1 and type 
2 AGN. Type 1 AGN show strong emission in the blue and UV 
regime of the spectrum, indicating a direct view to the accre- 
tion disk. Additionally, broad emission lines are present, which 
come from fast-moving material close to the central black 
hole. Type 2 objects do not show the blue/UV-bump. Emission 
lines are narrow, indicating that the material responsible for 
the emission is moving much slower and, thus, must be fur- 
ther away from the center. The difference between type 1 and 
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type 2 AGN is explained by obscuration of the central region. 
The obscuring structure is presumably an optically thick dusty 
torus, probably starting at the dust sublimation radius rsubi and 
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reaching out to 10 to 10^r,„bi (e.g., lGranato & Daneselll994t 
iBock et alJlzOOa |Radomski et alJ|2QQi. In type 1 AGN this 
torus is viewed approximately pole-on, while in type 2 objects 
it is seen close to edge-on. This implies that most type 2 objects 
harbor a type 1 nucleus - a suggestion which was confirmed by 
the observations of weak polarized broad emission lines hid- 
den by strong narrow lines in the Seyfert 2 galaxy NGC 1068 
('Miller & Antonucci'* 1983V The observed ratio of Seyfert 1 to 
Seyfert 2 galaxies of approximately 1 to 1.3 (Maiolino & Riek^ 
ll995l:lLacv et all2004l) suggests that the torus is both optically 
and geometrically thick. 

Additional evidence for the existence of a dust torus comes 
from the silicate feature at -lOyum in the spectral energy 
distribution (SED) of AGN. In type 1 AGN, this feature is 
expected to be detected in emission since the hot dust at 
the inner rim of the torus can be seen directly. Recent mid- 
infrared (MIR) spectra obtained with instruments on board the 
Spitzer satellite confirm the silicate emission feature in QSOs 
( Siebe nmorgen et al.l l2005t iHao et alJ l2005h and the LINER 
NGC 3998 (Sturm et al."2005^. For type 2 objects, the dust fea- 
ture is observed in absorption (e.g.. Jaffe et al. 2004) due to 
obscuration by the cold dust. Further evidence of a torus comes 
from near-infrared (NIR) bispectrum speckle interferometry 
(Iwittkowski et al.- 1998: Weigelt et al.^2004) and speckle inter- 
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ferometry d Weinberger et al.lfl999l) of the resolved nucl eus of 
NGC 1068. H- and 7^' -band images (IWeigelt et al.ll2004 show 
a resolved core structure of 1.3 x 2.8 pc. This resolved core 
is interpreted as dust close to the sublimation radius and an 
inner outflow cavity. Furthermore, strong evidence for a dust 
torus in NGC 1068 arises from long-baseline interferometry 
of NGC 10 68 with the V ery Large Telescope Interferometer 
(VLTI). .Ta fl^eetalJ (l2004 report on 8 - 13yum MIR observa- 
tions with the VLTI MIDI instrument. The wavelength depen- 
dence of the visibility inside and outside the silicate absorp- 
tion feature can be interpreted with a two-component model 
consisting of a hot (800 K) dust structure with an extension 
of < Ipc and a warm (320 K) structure of 3.4 x 2.1 pc. 
Finally, first near-infrared ^-band long-baseline interferometry 
of NGC 1068 wi th the VLTI and th e VINCI instrument was re- 
ported bv .Wittkowskiet al.l (l2004 . These VINCI observations 
detected a substructure with a size of < 3 mas (^ 0.2 pc). 

First pioneering radia tive transfer ni o deling of AG N dust 
tori was carried out by IPier & Krolil] (Il992[ h993h . They 
studied a torus with cylindrical geometry and a homogeneous 
smooth dust distribution. Further investigations with smooth 
dust distributions were reported bv lGranato & Danes3 (1 19941), 

t tathiou & Rowan-Robinson ( 1995) , an d iGranato et al.l 
97ft . Recently, Schartmann et al. (2005) presented very 
detailed radiative transfer modeling, introducing a physically 
mo tivated geometry and dust distribution (Camenzind 1995|). 

Krolik & BegelmanI |l988) discussed that smooth dust 
distributions probably cannot survive in the vicinity of an 
AGN. They proposed that the dust is arranged in clouds 
and that the observed velocity dispersion is actually a re- 
sult of their orbital motion within the gravitational poten- 
tial. First results of radiative tran sfer calculation s of su ch a 
clumpy medium were reported by iNenkova et al.l (l2002h and 
jDullemond & van Bemmel.(,2005|). They show that a torus con- 
sisting of clouds is capable of suppressing the strength of 
silicate emission and abs orption, as has been observed (e.g., 
[jaffe et al. 2004£ LSiebenmorgen et al. 2005: Hao et al. 2005: 
ISturm et al]l2005l) . It also became clear that the actual cloud 
distribution is a critical point. Therefore, physically motivated 
models of the torus are required. 

In this paper we present our 3D Monte Carlo radiative 
transfer modeling of clumpy AGN tori. We have developed a 
method which allows us to obtain high-resolution Monte Carlo 
simulations with relatively small computation times. We cal- 
culate the radiation fields of individual clouds at various dis- 
tances from the AGN and distribute these clouds within the 
torus region. After discussing the principle of our method, we 
present model spectra and images of clumpy tori and analyze 
the influence of various parameters. To illustrate the feasibility 
of our approach, we finally model the SED and visibilities of 
NGC 1068 simultaneously. 

2. Monte Carlo simulation of individual dust 
clouds 

This section gives an overview of the Monte Carlo radiative 
transf er treatment with our code mcsim_mpi ( Qh naka et alJ 
l2QQ5h and describes its application to dust clouds. In the fol- 




Fig. 1. Sketch of a clumpy torus. The hot and cold clouds are 
randomly distributed around a central accretion disk (blue). 
The size and optical thickness varies with radial distance. 

lowing, we present cloud SEDs and analyze some aspects of 
the Monte Carlo treatment. 

2.1. AGN primary radiation 

We consider a central super-massive black hole (smBH) which 
is fed by an accretion disk. The smBH and its accretion disk are 
surrounded by an optically and geometrically thick dust torus. 
JKrohk & Begelman (1988) argue that the temperature has to be 
-10^ K if the velocity despersion is of the order of 100 km s"^. 
This temperature greatly exceeds the dust sublimation tempera- 
ture. Thus, we arrange the dust in clouds. In Fig.[T] we illustrate 
the clumpy composition of the torus. 

The primary source of radiation is the accretion disk. We 
approximate the SED of the accretion disk by a broken power 
law: 

(A A< 0.03 yum 

AF;i oc I constant 0.03 jum < A<03 //m (1) 
[a-^ A> 0.3 jam 

It is illustrated in F ig. |2l The SED is a fit to QSO spectra 
(iManske et al.l 19981) with a characteristic steep decline towards 
the MIR. 

2.2. Geometry and dust cloud properties 

Our radiative transfer code mcsirn_mpi allows the model- 
ing of dust distributions of arbitrary 3 -dimensional geome- 
try. Since there is no knowledge about the shape of dust 
clouds around AGN, we model spherically symmetric clouds. 
For our calculations, we assume a standard galactic dust 
composition of 53% astronomical silicates and 47% graphite 
(Draine & Lee 1984) . The absorption and scattering efficien- 
cies Q^hs and Qsrp^ for MRN di stributions of grain sizes 
(|Mathis, Rumpl & Nordsieck' 1977) have been calculated with 
the DUSTY code (Ivezic et al. 1999) . In our Monte Carlo code, 
we use an averag e grain size of 0.1 /im representing the whole 
distribution (e.g. 'Efstathiou & Rowan -RobinsonI h994t IWolfl 
E003; Nenkovaet al. 2005). 

For the characterization of the cloud radiation field, three 
main parameters are needed: the bolometric AGN luminosity 
Lagn, the cloud's optical depth Tc/, and its distance r from the 



S. F. Honig et aL: 3D Radiative Transfer Modeling of Clumpy Tori and its Application to NGC 1068 



3 




0.01 0.10 1.00 

wavelength X []xm] 

Fig. 2. AGN source spectrum as used in our Monte Carlo sim- 
ulations. 



central source, r^/ is taken at the reference visual wavelength of 
0.55 ^m and denotes the optical depth through the cloud center. 
At other wavelengths, the optical depth is derived from gabs- 
Within a cloud, the dust is distributed homogeneously. 

2.3. Radiative Transfer Code 

Our radiative transfer code mcsim_mpi is based on the Monte 
Carlo approach by Bjorkman & Wood (200 11). A detai led de- 
scription of our code is given in Ohnaka et alJ Here 
we briefly summarize the concept. As for any Monte Carlo 
approach, the luminosity of the central radiation source is di- 
vided into monochromatic photon packages. If a photon pack- 
age hits a dust particle, it is either scattered or absorbed. In 
case the photon package is scattered, a random scattering an- 
gle is determined according to the difl'erential cross section. If 
it is absorbed, the dust particle is heated up and the photon is 
re-emitted at a frequency determined by its temperature. The 
approach by Bjorkman & Wood has the advantage that it is 
non-iterative for fixed opacities and thus requires less comput- 
ing time. As output, the code provides SEDs for scattered and 
absorbed photons, which are then added to a complete cloud 
SED. Furthermore, it is possible to obtain a number of cloud 
SEDs for diff'erent observer positions (=phase angles 0). 

The code has been tested intensively and was recently ap- 
plied to the modeling of VLTI MIDI visibilities of the evolved 
star IRAS 08002-3803 (Ohnak a et al. 2005). For testing our 
specific cloud geometry, we compared our cloud SEDs to re- 
sults obtained with the DUSTY code. DUSTY allows the mod- 
eling of a dust sphere in an external isotropic radiation field. 
Although this is not exactly the same geometry as our external 
illumination by a central source, we use the resulting SEDs for 
comparison with the hot-side SEDs of our mcsim_mpi clouds. 
While the overall SED shapes are quite similar (see Fig.O, the 
resulting diff'erences are expected due to diff'erent illumination 
patterns. 



100 



1 




X [\im] 

Fig. 3. Comparison of a hot-side cloud SED simulated with our 
Monte Carlo code mcsim_mpi (solid line) and a SED of a dust 
sphere with external isotropic illumination (obtained with the 
DUSTY code; dashed line). The Monte Carlo calculations were 
performed with 5 x 10^ photon packages (see Sect. 12.41) . 




Fig. 4. Temperature profile of a simulated cloud (r^/ = 150) 
with a high-resolution grid (150 x 32 cells; top) and a low- 
resolution grid (20 X 20 cells; bottom). The number of photon 
packages are 5x10^ for both simulations. The cloud was placed 
close to the dust sublimation radius, where the AGN is to the 
left. The hottest temperatures are ~1400K (white), the coolest 
are -200 K (dark red). 

2.4. Quality checks of Monte Carlo simulations 

Depending on Tc/, the temperature within a cloud decreases 
from the hot to the cold side. If the number of grid cells is 
too small, this temperature profile will not be reproduced prop- 
erly and parts of the cloud are either too cold or too hot. This 
specifically applies to clouds with high r^/ since the tempera- 
ture gradients can be steep. In Table [l] we show Monte Carlo 
temperature gradients for a cloud with r^/ = 100 at a distance 
ofr = Tsubi (T(rsnhi) = 1500 K) from the AGN. We fist the tem- 
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Fig. 5. SEDs of spherical dust clouds with optical depths of Td = 70 simulated with our Monte Carlo code rncsim_mpi. In the 
upper row, the panels represent clouds at distances r = rsubi and r = 2 rsubi; in the lower row they represent clouds at distances 
r = 4rsubi and r = 8 rsubi- For each cloud three SEDs are shown, as seen under phase angles of = 160° (green), = 110° 
(black), and cp = 60° (red). In terms of illumination, these phase angles correspond to 3%, 33%, and 75%, respectively. The 
absolute flux values are normalized to the cloud surfaces. 



Table 1. Temperature gradient variations for increasing num- 
bers of radial cells within a cloud (rhot = 1500 K). 



number 






radial cells 


[K] 


[K] 


10 


560 


660 


20 


580 


800 


40 


600 


980 


70 


610 


1090 


100 


640 


1160 


150 


640 


1180 



tion, the temperatures on the hot side are too cold compared to 
the high-resolution simulation. 

Furthermore, we have to ensure that enough photon pack- 
ages interact in each grid cell. As discussed above, we use 
100-150 grid cells in radial direction. Since the temperature 
gradients in the plane perpendicular to the radial direction do 
not need a very high resolution, we use only 32 cells. The 3- 
dimensional radiation field results from the rotational symme- 
try around the central radial axis. Thus, our clouds consist of 
~ 5x10^ grid cells. Tests with the code show that the results 
are quite good for 10^-10'^ photon packages interacting in the 
cloud grid cells. If we take less than ~ 10^ photon packages, 
the resulting SED becomes too noisy. 



perature gradient from the hot side to one tenth of the cloud di- 2 5 Cloud SEDs 
ameter inward, as well as the complete gradient from the front 

to the back side. As can be seen, the temperature gradients con- We performed calculations for a number of diff'erent cloud pa- 
verge once a certain number of radial cells is used. Therefore, rameters to study their SEDs and consequences for a clumpy 
it is quite important to use enough radial cells to sample the torus. For the simulations we used 10^ photon packages and 
temperature profile within a wide range of r^/ values. In Fig.|4] ~ 5 x 10^ grid cells for each cloud. 

we show the comparison of a cloud simulated with high and In Fig.|5]we show our results for clouds with r^/ = 70. The 

low resolution. For both simulations we used 5 x 10^ photon clouds were placed at diff'erent distances from the AGN (1,2, 

packages (see next paragraph). In the low-resolution simula- 4, and 8 rsubi). For each cloud we show SEDs for three diff'erent 
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phase angles 0. The phase angle is defined as the angle between 
the direction from the cloud to the observer and the direction 
from the cloud to the AGN, as seen from the cloud. = 0° 
corresponds to a cloud that appears fully illuminated, while 
= 180° represents the backside view of the cloud without 
direct AGN illumination. This concept is similar to the phases 
of the moon: depending on where the cloud is located around 
the AGN, only a fraction of its surface seen by the observer is 
directly illuminated by the AGN. 

The diff'erent panels in Fig. 13 show the dependence of the 
continuum peak on the distance to the AGN. The peak emis- 
sion of the SED continuum shifts to longer wavelengths for 
larger distances. In addition, the absolute flux decreases. At 
larger distances from the AGN, the relative flux of the silicate 
emission at -lOyum increases with respect to the continuum. 
Furthermore, it is remarkable that the silicate emission feature 
is not as pronounced as in torus model SEDs of smooth dust 
distributions. Although the SEDs of Fig 13 are model SEDs of 
individual clouds, the implications for a torus spectrum are evi- 
dent: because of the clumpy structure, the emission feature will 
not become very strong (see Sect.|4]and|3. 

It is interesting to follow the SED development from the hot 
frontside to the cold backside view of each cloud. The cloud 
with the smallest distance to the AGN (cf. upper left panel in 
Fig. 13 shows a clear shift of the continuum peak from around 
1.7 yum on the hot side to around 3 yum on the cold side. In ad- 
dition, the silicate feature in emission starts to decrease with 
respect to the continuum. This is a result of the temperature gra- 
dient within the cloud. On the hot side, the temperature is high 
enough to produce an emission feature. The emission basically 
takes place in the upper dust layers of the hot parts. When the 
cloud is seen from the cold side, the hot side's emission has 
to pass through colder layers within the cloud and sufl'ers from 
absorption. Furthermore, the colder layers produce less or no 
silicate emission. As a result, the overall feature decreases. At 
larger distances (cf. lower right panel in Fig. |3, the emission 
from the hot side barely compensates for the absorption within 
the cloud, so the silicate emission disappears in the SED of the 
cloud's cold side. In extreme cases of high r^/, large distances 
r, or no direct illumination, the emission may even turn into ab- 
sorption. This efl'ect is illustrated in Fig. |6l where SEDs from 
the cold side of clouds are shown for such cases. 



3. Torus assembly 

We will now present a simple parametrization of the torus. 
In the following, we describe our method for arranging the 
clouds within the torus. This method allows us to obtain model 
SEDs and model images of tori using the simulated cloud 
SEDs which have been calculated before. Afterwards, a phys- 
ical model for the cloud properties is outlined. To distinguish 
between properties of clouds and the torus, we will label all 
cloud parameters with the subscript cl, while non- sub scripted 
parameters always apply to the torus as a whole. 




1 10 
X [[Lm] 



Fig. 6. Model SEDs of the unilluminated side of dust clouds at 
two difl'erent distances from the central source. The silicate fea- 
ture is in absorption. The solid and the dashed SEDs represent 
clouds at 2 Tsubi with Td = 100 and Td = 400, respectively. The 
dashed-dotted SED comes from a cloud at 4 rsubi and Td = 250. 



3.1. Torus parameters 

To assemble the clumpy torus, it is necessary to get information 
about the cloud properties and the cloud distribution within the 
torus. A very simple approach to this problem is the model- 
ing of these properties with free parameters and correspond- 
ing power laws. Geometric torus input parameters are the total 
number of clouds within the torus A/^tot, their distribution r]r(r) 
and T](^((f) in radial and azimuthal direction, respectively, the 
scale height of the torus H(r), which also determines the cloud 
distribution in z-dire ction rj^ir, z), and its outer rim To (also see 
iNenkova etallbOOS!) . 

For the radial distribution, we use a power law T]r(r) oc 
r^. The z- and r-dependences of 7]^ are parameterized by a 
Gaussian-like ansatz: i.e., T],{r,z) = exp(z2/2H(r)2). This ac- 
counts for a higher cloud density towards the equatorial plane 
and flaring of the torus. For simplicity, r]^p((f) is set to unity so 
that 7/(r, z, (fi) = J]r(f) ' ^z(^, z) and J 7/(r, z = 0, ^) dr = 1. 

For convenience, we substitute the total number of clouds 
within the torus A/^tot by the number of cloud s A^o along the 
line o f view to the AGN in the equatorial plane (iNenkova et alJ 
I2OO5I) . The relation between A^tot and A^o is given by 

^ ^cl 

The basic torus parameters A^o. r](r,z,(f), H(r), and are 
supplemented by the properties of the clouds; in particular, 
their optical depth Td(r) and radius Rd(r)- Both parameters are 
written as r-dependent. Actually, whether Td changes with dis- 
tance or not is unclear. Theoretical arguments exist for a Td 
gradient (see Sect. 13.311 . Henceforth, a simple model with the 
power law parameters described above and listed in Table Sis 
called Sifree parameter model. 
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Table 2. Parameters of the free parameter model and the accretion scenario model. 



model parameter 


free parameter model 


accretion scenario 


number of clouds in equatorial LOV A^o 


free 


free 


optical depth Tci{r) 


oc r^^, Yj free 


oc M(r)^/V-3/2 


cloud size Rci{r) 


oc r^^, Ft? free 




cloud distribution 77^ (r) 


oc r^,a free 




scale height H{r) 


oc r^^, Th free 




outer radius Tq 


free 





3.2. Torus SED compilation and geometric 
peculiarities 

Once a distribution function for the clouds is selected, random 
positions for A/^tot clouds are determined according to T](r, z, (f). 
Every position is labeled with the corresponding cloud proper- 
ties Tci and Rci, and a template spectrum from the database is 
associated with this position. As it happens, some of the clouds 
are not directly illuminated by the AGN since another cloud 
is along the line of sight to the AGN. Such clouds, however, 
can be indirectly illuminated by the radiation coming from sur- 
rounding clouds. Therefore, we must distinguish between these 
indirectly illuminated clouds (Second Order Clouds; SOCs) 
and clouds that are directly illuminated by the AGN (First 
Order Clouds; FOCs). 

The treatment of FOCs and SOCs is quite different. For 
FOCs we simply take the SEDs from the Monte Carlo output, 
as described in Sect. l2.5l Since there is no direct AGN illumina- 
tion for a SOC, one has to approximate the illumination caused 
by all other clouds. Given the grid resolution that is needed for 
proper handling of the clumps, the necessary distances between 
the clouds, and the total number of clouds of several 10"^, an ap- 
proximation cannot be avoided since, otherwise, the number of 
photon packages would be 10^^ - 10^^ and more (see Sect. l2.4li . 
This is difficult to achieve with current computation power. 

In general, a SOC is located within the radiation field of 
the surrounding clouds. To approx imate this diffu s e radi ation, 
we follow the scheme proposed by Nenkova et "Zl fcOO.Sh . The 
main contribution to a SOCs ambient radiation comes from the 
hot side of FOCs in the vincinity of the SOC. If we consider the 
diffuse radiation field to be isotropic, it is possible to approxi- 
mate its SED at a given distance r from the AGN by averaging 
a FOC SED over all phase angles 0. To finally obtain the SED 
of a SOC, the average FOC spectrum can be taken as the in- 
put spectrum for an ambient radiation field in our Monte Carlo 
code. The stength of the ambient radiation field is derived by 
taking into account the fraction of the sky covered by FOCs as 
seen by a SOC near the sublimation radius. Because the am- 
bient radiation field is assumed to be isotropic, the simulated 
SOC shows an identical SED for any phase angle 0. The results 
obtained from these calculations form the database for SOC 
templates. The method is discussed in AppendixlAlin more de- 
tails. 

When all clouds have been associated with template SEDs 
from the FOC and SOC database, the torus SED can be com- 
piled. For this purpose, a torus inclination angle / has to be 
chosen. A torus seen pole-on has / = 0°, while the edge-on 



view corresponds to / = 90°. As predicted by the Unification 
Scheme, type 1 AGN should have small /, whereas / should be 
larger (e.g., > 45°) for type 2 objects. The inclination / also 
defines a line of view (LOV) from the observer to the torus 
clouds. With respect to this LOV, each cloud is seen under a 
specific phase angle 0, depending on its position in the torus. 
The corresponding SED as seen by the observer is taken from 
the templates, and the LOV of each cloud is checked for obscur- 
ing clouds along the LOV. If there are other clouds, the tem- 
plate SED is attenuated according to the wavelength-dependent 
Tc/ -profile and the fraction of coverage. The resulting SED is 
added to the total torus SED. Once all clouds are sifted through, 
the torus SED is complete. Since this method of compiling 
the torus spectrum does not demand any special symmetry, it 
is also applicable to non-axisymmetric cloud distributions or 
other geometries. Actually, the statistical process of distribut- 
ing clouds already introduces 3 -dimensional asymmetries. 

3.3. Physical torus models and cloud properties 

Power law models, as described in Sect. l3.ll are widely used in 
torus studies, although it would be desirable to have the param- 
eters fi xed by theoretical arguments. For example. ICamenzindl 
(Il995h describes a smooth dust torus where the dust is pro- 
duced in stars of a circumnuclear stellar cluster (cnSC). The 
toroidal geometry is a result of the effectiv e gravitational po- 
tential coming from the smBH and the cnSC (ISchartmann et alJ 
I2OO5). 

A ph ysically motiva t ed mo del f or clumpy dust tori was pro - 
posed bvVoll mer et alJ (EoO?) and'Becke rt & Duschll (I2QQ4 . 
Its basic idea is that the dust clouds within the torus are self- 
gravitating. Their size is close to the shear limit given by the 
gravitational potential of the cnSC and the smBH. Actually, the 
torus is in a region where the influences of both the cnSC and 
the smBH are important. The cnSC is supposed to be spher- 
ically symmetric and isothermal; i.e., the relation of the en- 
closed mass of the cnSC is given by M^^(r) oc r. In this torus 
model, the material that forms the torus is accreted from large 
distances of more than lOOrsubi- The accretion scenario as- 
sumes effective angular momentum transfer due to cloud colli- 
sions. This leads to large scale heights for high accretion rates 
above the Eddington limit of the central smBH. Since the en- 
closed mass in the torus is a combination of smBH and cnSC 
mass, the Eddington limit within the torus is higher than in the 
accretion disk. An advantage of this model is that it reduces the 
number of free torus parameters. In addition to the AGN bolo- 
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Fig. 7. Torus SEDs of 10 different random cloud arrangements for the set of parameters listed in Table [S] The panels in the first 
row show SEDs for inclinations of / = 90° (edge-on; left) and / = 45° (right); the panels in the second row show / = 30° (left) 
and / = 0° (pole-on; right). 



metric luminosity Lagn and its black hole mass Mbh, one has 
to specify the mass accretion rate M through the torus and the 
properties of the cnSC. The latter is defined by its core radius 
and the core mass M^^^. The dependencies of the accre- 
tion scenario are listed in Table |2l in comparison to the free 
parameter model and constraints coming from observations. In 
the table, M(r) denotes the sum of cnSC and smBH masses. 
Furthermore, the accretion scenario doesn't need an outer ra- 
dius for the torus. Since it predicts Rd oc r^l^ • M(r)"^/^, and 
Mci oc r^l^ • M(r)"^/^, the optical depth of a cloud is given by 

^cl 

For a n isothermal cnSC, a s observed in the center of our galaxy 
(e.g. ISchodel et alJ 12003). M(r) is proportional to r so that 
Tci(r) oc r"^ This implies that Tc/(r) steadily decreases with in- 
creasing distance from the AGN. At larger distances, the clouds 
become larger and less dense; they smoothly pass into the sur- 
rounding galactic ISM. 

Some interesting predictions of the accretion scenario are 
the cloud sizes and masses. The total mass for clouds - gas and 
dust - close to the sublimation radius rsubi should be of the or- 
der of 1 ... 5 M©, while the cloud radius is approximately 0.01 
to 0.05 pc. Since Rd and Md weakly depend on the enclosed 
mass, these values are similar for a wide range of smBH and 



Table 3. Parameters as used in our parameter studies. 



Torus 


Model Value 


Clouds 


Model Value 


^AGN 
^subl 


2x Li2 
0.8 pc 


Rciir) 
No 


0.02pcx(^) 

^ V '•subl / 

8 

/ \-0.8 
40x(T7r^) 

VlOrsubl/ 


d 

rjrir) oc 


lOMpc 
a = 1.5 


Tcl 


H(r) 




dust 


see Sect.lZ2l 


ro 


70 Tsubl 







cnSC masses. Henceforth, for the modeling we use cloud prop- 
erties which are consistent with the accretion scenario. 



4. Results 

In this section we present model SEDs and model images ob- 
tained from our torus simulations. The eff'ect of the clumpy 
structure on the torus is analyzed by parameter studies. We 
briefly discuss the impact of dumpiness on the modeling of 
observations. For our parameter studies, we have chosen pa- 
rameters listed in Table [3l Changes and variations to this set of 
parameters are noted in the text. 
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Fig. 8. L-band model images of the torus described by the parameters listed in Table |3l (same parameters as used for the SED 
calculations shown in Fig.Q). Left column: images obtained by averaging model images of -200 different random cloud arrange- 
ments. Middle and right columns: model images for two particular random cloud arrangements. From top to bottom: / = 90°, 
45°, 30°, and 0°. 
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Fig. 10. The dependence of the SED on rir (model parameters in Table [3. The upper panels show edge-on SEDs for 10 different 
random cloud arrangements for a = 1.0 (left) and a = 1.5 (right). In the lower left panel we present the edge-on view for a = 2.0, 
while the lower right panel combines pole-on views of the three models with a = l.O (top), a = 1.5 (middle), and a = 2.0 
(bottom). 




Fig. 9. ^-band model image for non-isotropic, sinusoidal ra- 
diation from the accretion disk for one particular random cloud 
arrangement, edge-on view, and model parameters listed in 
TableEl 

4.1. Statistical cloud distribution 

The distribution of clouds follows a statistical process. Thus, it 
is expected that different statistical realizations of this distribu- 
tion will change the torus SED. Since recent studies of clumpy 



tori use statistical averages (iNenkova et al.ll2002L Eool . It IS 
also worth investigating the scatter that has to be expected. In 
Fig. Q we show the SED variation for one set of parameters 
(listed in Table 13. The only difference between the SEDs is 
the random arrangement of the clouds, following the same dis- 
tribution function 7/(r,z, (fi). Analyzing the data shows that the 
slope of the NIR flux is comparable for all random cloud ar- 
rangements. The absolute fluxes and the silicate feature differ 
by a factor of 2 - 3 for the edge-on view. The large SED scatter 
at / = 30° is a result of the transition from edge-on to pole-on 
view. The FIR fluxes get close to each other at longer wave- 
lengths since the optical depth of the torus becomes smaller 
than unity, so that obscuration no longer plays an important 
role. 

We investigated the origin of the dependence of the SED on 
the random cloud arrangement. As we found out, these varia- 
tions are caused by the FOCs. They are the main contributors at 
the shorter wavelengths of the SED. A few bright FOCs at high 
latitudes can change the SED considerably since high-latitude 
clouds are usually less obscured. 

In Fig. [S] we present model images for inclinations / = 90°, 
45°, 30°, and 0°. For each inclination angle, we show two 
different kinds of images: the left column presents an aver- 
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Fig. 11. SEDs for different flaring (see text; model parameters in Table [Sj. From upper left to lower right: H oc r^-^, r^-^, r^-^, 
and r^-^, respectively. Each panel consists of SEDs for 10 difl'erent random cloud arrangements seen edge-on (black) and pole-on 
(red). 



age image derived from -200 difl'erent random cloud arrange- 
ments; the mid and right columns show images calculated for 
two individual random cloud arrangements (wavelength range 
4.65 -4.95 yum; for model parameters as listed in Table O. In 
all images, the assumed faint emission from the central AGN 
is not shown (see Sect. 12. II and Fig.|2|). The images obtained 
for particular random cloud arrangements show a similar over- 
all appearance with distinct diff'erences in the substructure. The 
grainy structure in the images is a result of the dumpiness of 
the torus and not related to noise of the radiative transfer simu- 
lations. It illustrates the appearance of a clumpy torus for very 
high-angular resolution observations. 

In the average image for / = 90°, a polar cavity with an X- 
shaped structure is visible. At / = 45° and 90°, a bright inner 
structure is the dominant source of emission, which is caused 
by the inner torus rim. At / = 30°, parts of the inner rim of the 
torus can be seen as a crescent structure. The central dark re- 
gion for / = 30° and 0° is the inner dust- free region. The appar- 
ently well-defined brightness steps in the averaged images will 
disappear as more diff'erent clouds are added to the database of 
SOCs. 

The model images presented in Fig. [S] were calculated for 
an isotropic radiation of the accretion disk. Fig. |9| shows a 
model image obtained for a non-isotropic, sinusoidal radiation 



field of the accretion disk. If the torus is illuminated by this 
non-isotropic radiation, a dark dust lane appears in the disk 
plane when the torus is viewed edge-on. 

4.2. Radial Distribution function 

To illustrate the eff'ect of the radial cloud distribution function 
T]r(r) on the SED, we keep all parameters in Table Inconstant, 
except for exponent a of the radial power law. Fig.[lO|shows the 
resulting SEDs for a = 1.0, a = 1.5, and a = 2.0 for both the 
edge-on and pole-on views of the torus. The relative depth of 
the silicate feature in absorption becomes larger with growing 
a. Additionally, the NIR flux decreases for both the pole-on and 
edge-on views, and the variations become stronger. This can be 
explained by the higher density of clouds with large r^/ close 
to Tsubi- As a result, the obscuration of the FOCs in that area 
becomes higher and their total flux decreases. Small difl'erences 
in the arrangement may then lead to higher variations. 

It is interesting to look at the evolution of the silicate emis- 
sion feature in the pole-on view. The feature is quite flat and 
less pronounced than for smooth dust distributions, and it de- 
creases for larger a. This efl'ect comes from the higher cloud 
density of FOCs close to rsubi when emission and absorption in 
that region become equally important. Tori with a^2 can flat- 
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Fig. 12. SEDs for outer radii of - 70 rsubi (black) and To = 
30 Tsubi (red) seen edge-on and pole-on (model parameters in 
TableI3. 

ten out the emission feature for small / or even show moderate 
absorption. 

4.3. Jhe effect of flaring 

The z-dependence of the cloud distribution function is a 
Gaussian profile rj^ ^ exp (-z^/2//(r)^) (see Sect. 13. 111 . Thus, 
the toroidal structure flares depending on the r-dependence of 
H. In Fig. [n] we show the model SEDs for H oc r^-^, r^-^, r^-^, 
and r^-^. For the edge-on geometry, the flaring has some efl'ect 
on the relative depth of the silicate feature. With small flaring 
the feature is deeper, and the overall flux is also lower by a fac- 
tor of ~3 - 4. This can be explained by the absence of clouds at 
higher latitudes which are usually less obscured and thus con- 
tribute significantly to the NIR and MIR fluxes. 

For the pole-on view, the importance of flaring is more evi- 
dent. The NIR and MIR fluxes decrease significantly when the 
flaring becomes stronger since there are more and more clouds 
at high latitudes that obscure the denser parts of the torus struc- 
ture. As a result, the emission feature vanishes for H oc r^-^ and 
can even turn into weak absorption for stronger flaring. 

4.4. The outer radius of the torus 

Since the main torus emission in the NIR and MIR comes from 
the innermost part, it is expected that the SED is less sensi- 
tive to the outer boundary of the torus. This especially applies 
if Tci decreases with r, as proposed by the accretion scenario 
(Sect. 13.3b . In this case, the innermost part of the torus dom- 
inates since absorption by more distant clouds is weak. Thus, 
the strongest efl'ect of the outer torus radius is expected in 
cases where distant clouds significantly contribute to the ob- 
scuration of the inner torus. To test this, we investigated the 
influence of Vq for the extreme case of r^/ = 70 (constant) 
throughout the torus. We simulated SEDs for Vq = 30 rsubi and 
To = 70 Tsubb respectively. 

Our results are shown in Fig. [El The overall shape of the 
SED is similar for both model To, while there are some dif- 



Table 4. Photometric data for the nuclear region of NGC 1068 
from high-resolution observations. 



Wavelength 


Flux 


Phot. Aperture 


Reference 


fim 


Jy 


arcsec 




1.3 


0.00048 


0.2 


R98 


1.6 


0.07 + 0.02 


0.1 


W04 


1.8 


0.0054 


0.2 


R98 


2.1 


0.35 + 0.09 


0.1 


W04 


2.2 


0.056 


0.27 


P04 


3.5 


1.34 


0.4 


MOO 


3.8 


1.42 


0.4 


M03 


4.5 


2.5 


0.27 


P04 


4.7 


2.72 + 0.14 


0.4 


M03 




j.Zj ± U.jZ 




A/rnn 

iVlUU 


7.7 


8.8 + 1.0 


0.4 


TOl 


8.7 


10.0+1.2 


0.4 


TOl 


9.7 


7.04 + 0.80 


0.4 


TOl 


10.4 


7.94 + 0.87 


0.4 


TOl 


11.7 


17.8 + 2.2 


0.4 


TOl 


12.3 


17.2 + 2.1 


0.4 


TOl 


18.5 


20.2 + 3.4 


0.4 


TOl 



R98 fR ouan et al. 1998 ). W04 fWeigel t et al. 20041 



P04 (Prieto, published in Schartmann et al. 2005), 

MOO rMarco & Alloin 2000) . M03 rMarco & Brooks 2003) . 

TOl rTomonoetal.2001) 

ferences in the absolute fluxes. Since the central region is less 
obscured for To = 30 rsubi, the NIR and MIR fluxes are higher. 
On the other hand, the FIR flux of the To = 70 rgubi model is 
above the flux for the smaller radius since the outer cold part is 
missing in the latter model. 

5. Simultaneous modeling of the SED and 
visibilities of NGC1068 

The goal of our radiative transfer modeling studies of 
NGC 1068 is to find a clumpy torus model which can simulta- 
neously reproduce both (la,b) the observed high spatial resolu- 
tion IR SED and (2a-c) the interferometric visibilities obtained 
during the last few years: 

(la) SED values in the wavelength range of 1.3 to 18 jim 
Hsted in Tablegl 

(lb) high angular resolution MID IR SED in the wave- 
length range from 8 to 13 urn obtained using the VLTI MIDI 
instrument ( Jaff'e et al. 2004), 

(2a) speckle int erferometric H- and /sT-band visibilities 
(|Weigeltetal.||2QQi, 

(2b) MIR visi bilities (8 - 13/i m) obtained with the VLTI 
MIDI instrument llafl^e et aJ2004h . 

(2c) i^-band visibilities obtained with the VLTI VINCI in- 
strument ( Wittkowskiet al. 2004). 

We tried to find a set of model parameters that can simul- 
taneously reproduce all of these observations. For this goal we 
scanned the parameter space, as described in Sect. |4j and ob- 
tained the clumpy torus model parameters listed in Table 13 
Figs. fT3l and fT6l present the model SED and M-, and N- 
band model images. 
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Fig. 13. Comparison between the observed high-resolution 
SED of NGC 1068 (blue diamonds and red curve) and our 
model SED (model param eters in Table 5). The red line rep- 
resents the total MIDI flux ('.Ta fl^eetal.lEool . The shaded area 
shows the range of model SED variations obtained for 10 dif- 
ferent random cloud arrangements. The dark grey line is the 
average of the 10 individual model SEDs. 

Table 5. Parameters of our NGC 1068 clumpy torus model. 




Baseline [m] 

Fig. 14. Comparison between the ^-band model visibilities 
(grey lines; model parameters of Table ISj for 10 difl'erent ran- 
dom cloud arrangements and the observed visibilities (red sym- 
bols) at 0-6 m and 46 m baseline at PA 45°. The green dashed- 
dotted line shows a Gaussian profile with 28 mas FWHM (=2.0 
pc; see Table |6|). 

Table 6. Comparison between the observed interferometric 
sizes and clumpy torus model sizes. 



Torus 


Model Value 


Clouds 


^AGN 


0.24 X Li2 


Rciir) 


^subl 


0.28 pc 


No 


d 


14.4 Mpc 


Tcl 


rjrir) oc 


a = 1.5 




H(r) 




dust 


ro 


70 Tsubl 




i 


55° 





O.Olpcx(^) 
10 

V 10 rsubl / 

see Sect. 12.21 



TT" 



Fig.fT3l shows a comparison of NGC 1068's average model 
SED, the variation range of the model SEDs obtained for 10 
diff'erent random cloud arrangements, and the observed high 
angular resolution SED. There is a satisfactory agreement be- 
tween the observed and modeled SED. 

Fig. [21 shows the visibilities of our NGC 1068 model 
(Table 13) in the range of 8 - 13 yum. The model visibilities are 
roughly in agreement wit h the observed M IDI visibilities for 
42 m and 78 m baselines ( Jaff'e et al. '2004'). The discrepancies 
are possibly caused by the diff'erence between our assumed dust 
properties (see Sect. 12.21) and the real ones. 

In addition to the speckle and MIDI observations, it was 
possible to observe NGC 1068 with the V LTI VINCI instru- 
ment in the ^-band (Iwittkowski et all2004l) . The obtained vis- 
ibility of -0.4 at a projected baseline of 46 m corresponds to a 
structure that is smaller than 3 mas or ^0.2 pc; i.e., at least ~ 10 
times smaller than the resolved ~30 mas core. The large visi- 
bility of 0.4 at 46 m baseline can be explained by at least two 
diff'erent types of structures within the 18 x 39 mas core: (1) a 
single compact object (in addition to the 18 x 39 mas core) with 
a diameter in the range of and 3 mas, or (2) several objects 
with diameters in the range of and 3 mas. Since there is only 



Band 


model size 


interf. size 


Reference 


H 


2.1+0.5 pc 


1.3 X 3.1 pc 


W04 


K 


2.0+0.3 pc 


1.3 X 2.7 pc 


W04 


N 


2.7+1.0 pc 


<2.1 x3.4pci 


J04 



W04 fWeiselt et al. 2004), J04 Qaffe et al. 2004) 

^ warm component: 2.1 X 3.4 pc; hot component: < 1 pc 



one visibility point currently available in the K band at long 
baselines, it is not yet possible to decide whether this substruc- 
ture consists of one or several ^0.2 pc objects. One possible ex- 
planation for this substructure is graininess due to dumpiness, 
as seen in the images in Figs. [S] and [l6l The corresponding K- 
band visibilities are shown in Fig.[^ Another possible expla- 
nation is that we see the central ac cretion flow viewed thr ough 
only moderate extinction directly (IWittkowski et al.l2004l) . 

Fig. [161 presents model images of our clumpy torus model 
of NGC 1068 (see Table 13 . These model images are calcu- 
lated for the wavelength ranges 2.1 - 2.3 yum, 4.7 - 5.0yum, 
and 9.5-10.3 yum, which approximately correspond to 
and A/^-bands, respectively. The images in the upper row are ob- 
tained for one particular random cloud arrangement, whereas 
the images in the lower row are averages over ~200 difl'erent 
random cloud arrangements. The grainy structure in the im- 
ages in the upper row is caused by the dumpiness of the torus 
and is characteristic for the clumpy torus structure. 

Finally, it is possible to compare the resolved //-, and 
N-band sizes of the compact core of NGC 1068 with the size 
predictions of our model described in Table 13 From the model 
images, we derived FWHM diameters (Gaussian model fits) 
of 2.1, 2.0, and 2.7 pc for the and A^-band model di- 

ameters. The comparison (Table 13 with the interf erometri- 
cally observed //-, and A/^-band diameters (1.3 x 3.1 pc. 
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Fig. 15. Comparison between the observed NGC 1068 MIDI visibilities (red line) at projected baselines of 42 m (upper panel) 
and 78 m (lower panel), and our model visibilities (model parameters in Table [Sj. The grey lines show the visibilities for 10 
different random cloud arrangements. The blue line represents the average model visibility. 



1.3 X 2.8 pc, and < 2.1 x 3.4 pc, respectively) shows that the 
observed and modeled diameters are roughly in agreement. In 
other words, our model allows us to simultaneously reproduce 
both the observed SED and the and A/^-band diameters 

of NGC 1068. 

6. Summary 

We presented a new two-step method for 3 -dimensional radia- 
tive transfer modeling of clumpy dust tori of AGN. First, SEDs 
of individual dust clouds are simulated with our Monte Carlo 
code. Then, the modeled clouds are randomly arranged within 
the torus region. The cloud properties are derive d from a phys- 
ically motivated accretion scenario proposed bv^ VoUmer et al.l 
(l2004h and Beckert & Duschl ( 2004 ). We investigated the influ- 
ence of clumpy dust distributions on model SEDs and images 
and presented model calculations for NGC 1068. 

For type 2 objects, we find a moderately pronounced ab- 
sorption feature which is in agreement with observations. Type 
1 model SEDs show moderate emission. The modeled silicate 
emission-to-continuum flux ratio of ~ 2 (see Figs. [El and 
fTTl agrees with recently observed emissi on features in QSQs 
and a LINER where the ration is 1.5-2.5 (ISieben morgen et al.' 
a005: Hao et al. 2005; Sturm et al. 2005|). In addition, the to- 
tal NIR and MIR model fluxes are approximately in agreement 
with observations. It is not necessary to assume additional dust 
in the polar region to reproduce the observed fluxes. 

The actual random cloud arrangement has a significant ef- 
fect on the SED and images. For the same set of large-scale 
torus parameters, we find remarkably diff'erent SEDs for diff'er- 



ent random cloud arrangements. The depth of the silicate fea- 
ture decreases if the number of unobscured clouds is increased 
in comparison to other random arrangements. Strong statistical 
variations of the optical depth of the torus due to the dumpiness 
have already been predicted bv lNenkova et al.l (|2002). They ar- 
gue that due to Poisson statistics, dumpiness is able to account 
for the high variation in X-ray column densities in Seyfert 2 
galaxies. We show that dumpiness not only afl'ects absorption 
of radiation from the central source but also has strong conse- 
quences for the IR emission of the torus itself. 

For type 1 AGN, we find that the silicate emission may van- 
ish if the toroidal geometry has strong flaring or a high cloud 
density towards the center. This may apply to Seyfert 1 nuclei 
with no emission feature observed. In such cases, the absence 
of a silicate feature does not necessarily mean that there is no 
dust torus present. 

In Sect. 13 we presented our modeling of the nucleus of 
NGC 1068. We showed that our clumpy torus model allows us 
to roughly reproduce both NGC 1068's SED and visibilities in 
the NIR and MIR simultaneously. The obtained strong graini- 
ness in the images is caused by the dumpiness. We also showed 
that clumpy torus models are able to reproduce deeper silicate 
absorption features with growing interferometric baseline, as 
reported bv Llafl^e et alJ too4 . 
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Fig. 16. NGC 1068 model images of our clumpy torus model described in Table|5l From left to right: M-, and A/^-band images 
for an inclination angle of / = 55°. The upper row shows images obtained for one particular random cloud arrangement, whereas 
the images in the lower row are averages of over -200 different random cloud arrangements. The grainy structure in the images 
in the upper row is caused by the dumpiness of the torus. 



Appendix A: Discussion and assessments of the 
method 

As described in l3.2l our treatment of the SOC heating assumes 
that we have the same number of FOCs around a SOC at all dis- 
tances from the AGN. In reality, the number of FOCs, however, 
decreases with growing distance. As a result, we overestimate 
the ambient radiation field that heats the SOCs. We discuss the 
influence of our treatment of the ambient radiation field on the 
torus simulations in the following sections. 

A.1. Energy conservation 

One important test of our simulations is the study of energy 
conservation. It is possible to derive a theoretically expected 
luminosity of the torus. In our model, the average number of 
clouds along a LOS, (A^), is given by a Gaussian distribution: 

(A.) = ;Vo.exp(-i^) = ;Vo.exp(-l^) (A.l) 

where a is the latitudal angle in spherical coordinates (z = r • 
tan Of) and H = H{r) is the scale height. For simplicity, we 
assume H{r) = • ^ in the following discussion. (A^) is the 



expectation value for Poisson statistics (iNenkova et al.ll200'2l) . 
The probability of finding n clouds along a LOS is given by 

(NT 

P(^)= ^exp(-<A^)). (A.2) 
n\ 

The probability that we see the AGN unobscured along a LOS 
is obtained for ^=0: 

P(0) = exp(-<A^)) = exp -A^o • ^^P -2^;2 ^^-^^ 

It is our aim to find the probability that a photon emitted by the 
AGN in a random direction is hitting a cloud in the torus. For 
this goal, we simply determine from eq. (IA.3II the probability 
Wo that a photon emitted by the AGN in a random direction is 
escaping the AGN undisturbed by averaging all possible LOS: 

2 



The probability wt that a photon hits a cloud in the torus is 
given by 

Wt = I - Wo (A.5) 

Thus, the theoretically expected luminosity Ltheor of the torus 
is 

^theor = • LagN = (1 " >^o) ' ^AGN (A.6) 
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This is the theoretical reference luminosity for our simulations. 
For energy conservation in the Monte Carlo simulations of the 
clouds, Ltheor should be equal to the (unobscured) re-emitted 
luminosity Lpoc of all first order clouds in the torus. We made 
a number of simulations which showed that in our simulations 
^theor = ^Foc to an accuracy of ~2%. 

Finally, this luminosity Lpoc should be re-emitted by our 
simulated tori when obscuration and SOCs (ambient radiation 
field) are considered in the simulations, so that Ltheor = Lpoc = 
Lsim- Our studies showed that in our simulations 

Ltheor = Lf0C± ~ 2% = Lsim± ~ 5%. 

Differences of the order of -5% are observed, depending on 
the actual cloud distribution function. From this analysis we 
conclude that energy is conserved in our simulations with sat- 
isfying accuracy. 

k.2. SED shape 

We have carefully investigated the eff'ect of our SOC approx- 
imation on the total torus SED shape. Our SOCs are gener- 
ally too hot since the heating by surrounding clouds is over- 
estimated (see Sect. 13.21) . The ambient radiation field around a 
SOC in the torus is determined by considering the fraction of 
the sky around a SOC that is covered by FOCs. This fraction is 
derived within the inner 2.5 sublimation radii, where the heat- 
ing of SOCs by surrounding FOCs is strongest within the torus. 
This treatment results in an overestimation of the ambient radi- 
ation field around SOCs at larger distances from the AGN. 

We investigated at what wavelengths our SOCs contribute 
to the torus SED. By removing the SOC emission from our 
simulations, we found out that the SOC contribution is rele- 
vant only longward of ~ 20yum. In the NIR and MIR up to 
20 yum, there is no diff'erence in the torus SED simulations with 
and without SOCs. The discrepancy between a torus simulation 
with SOCs and without SOCs becomes a factor of ~2.5 only 
longward of lOOyum. This factor is constant for wavelengths 
larger than 100 jim. 

In a next step, we investigated the influence of the FOCs 
on the torus SED. In our current treatment, we only consider 
direct heating by the AGN. FOCs, however, are also heated by 
surrounding clouds. This missing heating of FOCs mainly con- 
cerns the cold backside of the FOCs. To get an estimation of 
this additional heating, we added a SOC SED from our SOC 
treatment (i.e., an overestimated SOC SED) to the FOC back- 
side SED. We found out that in our simulations the contribution 
from additionally heated FOC backsides is small at all wave- 
lengths. In the NIR and MIR, no difl'erence in the SED is seen 
at all. In the FIR, a small increase in flux can be observed. The 
increase at lOOyum is -10%. Thus, the contribution of heated 
FOC backsides can be considered as insignificant in our simu- 
lations. 

The SEDs obtained with our SOC treatment represent an 
upper limit for the real SED. In addition, our simulations with- 
out SOCs set a lower limit for the SED. The real SED is some- 
where in-between. From the discussion above, we can set upper 
limits for the maximum overestimation and conclude: 



- In the short wavelength range up to 20 yum, our SOC treat- 
ment has almost no influence on the SED (i.e., error <2%). 

- Between 20 yum and 30 yum, our SOC treatment can cause 
maximum deviations from the real SED of the same order 
as the variations by difl'erent random cloud arrangements 
(i.e, -5-10%). 

- Longward of lOOyum, the flux overestimation due to our 
SOC treatment is less than a factor of ~2.5. 

- The influence of our SOC treatment on the luminosity is 
less than 5% (see Sect. lAJTi . 

- The heating of FOC backsides by surrounding clouds is in- 
significant for NIR and MIR, and is only -10% at 100 yum. 

The studies also highlight that the hot front sides of the FOCs 
are the dominant contributors to the ambient radiation field that 
heats the SOCs. Therefore, we conclude that our simulations 
are satisfying shortward of 20 yum. An improvement of the am- 
bient radiation field for longer wavelengths by simulating the 
local FOC-to-SOC-relation and the r-dependence of the FOC 
density for any cloud distribution will be the topic of a future 
paper. 
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